Phylogeography using mitogenomes: A rare Dipodidae, Sicista betulina, in North‐western Europe

Abstract Repeated climatic and vegetation changes during the Pleistocene have shaped biodiversity in Northern Europe including Denmark. The Northern Birch Mouse (Sicista betulina) was one of the first small rodent species to colonize Denmark after the Late Glacial Maximum. This study analyses complete mitochondrial genomes and two nuclear genes of the Northern Birch Mouse to investigate the phylogeographical pattern in North‐western Europe and test whether the species colonized Denmark through several colonization events. The latter was prompt by (i) the present‐day distinct northern and southern Danish distribution and (ii) the subfossil record of Northern Birch Mouse, supporting early Weichselian colonization. Samples from Denmark, Norway, Sweden, Russia, Latvia, Estonia, and Slovakia were included. Mitogenomes were obtained from 54 individuals, all representing unique mitogenomes supporting high genetic variation. Bayesian phylogenetic analysis identified two distinct evolutionary linages in Northern Europe diverging within the Elster glaciation period. The results of the two nuclear genomes showed lower genetic differentiation but supported the same evolutionary history. This suggests an allopatric origin of the clades followed by secondary contact. Individuals from southern Denmark were only found in one clade, while individuals from other areas, including northern Denmark, were represented in both clades. Nevertheless, we found no evidence for repeated colonization's explaining the observed fragmented distribution of the species today. The results indicated that the mitogenome pattern of the Northern Birch Mouse population in southern Denmark was either (i) due to the population being founded from northern Denmark, (ii) a result of climatic and anthropogenic effects reducing population size increasing genetic drift or (iii) caused by sampling bias.

were obtained from 54 individuals, all representing unique mitogenomes supporting high genetic variation. Bayesian phylogenetic analysis identified two distinct evolutionary linages in Northern Europe diverging within the Elster glaciation period. The results of the two nuclear genomes showed lower genetic differentiation but supported the same evolutionary history. This suggests an allopatric origin of the clades followed by secondary contact. Individuals from southern Denmark were only found in one clade, while individuals from other areas, including northern Denmark, were represented in both clades. Nevertheless, we found no evidence for repeated colonization's explaining the observed fragmented distribution of the species today. The results indicated that the mitogenome pattern of the Northern Birch Mouse population in southern Denmark was either (i) due to the population being founded from northern Denmark, (ii) a result of climatic and anthropogenic effects reducing population size increasing genetic drift or (iii) caused by sampling bias.

K E Y W O R D S
divergence time, genetic diversity, mitogenomes, Phylogeography, population structure, Sicista betulina

| INTRODUC TI ON
The present species pool and biogeography of European plant and animal populations are strongly shaped by repeated climate changes during the Pleistocene epoch and by farming practices during the mid-and late Holocene (Emanuelsson, 2009;Lang, 1994). These climate changes have impoverished the rich Miocene and Pliocene European flora (Mai, 1995;Svenning, 2003) and over much of the northern Europe extensive Pleistocene ice sheets (Ehlers & Gibbard, 2004) have resulted in recurrent "tabula rasa" situations.
Until the Latest Glacial Maximum (LGM) around 20,000 years ago, northern Europe including parts of Denmark was covered by ice. During the Late-and Post-glacial periods when the glaciers retreated, a succession of habitat types from tundra to deciduous and coniferous forests occurred. Immigration waves of colonizing fauna started as early as in the last part of the de-glaciation stage (c. 17,000-15,000 years ago; Figure 1a). In the beginning, immigration of faunas from east and south was facilitated by a continuous continent; however, due to rising sea water levels and lowered landmasses, around 8000 years ago, continental Denmark was transformed into a group of islands and one peninsula (Aaris-Sørensen, 2016), creating barriers to further immigration leading to local extinctions. Britain became an island, as the land between Denmark and Britain (Doggerland) was flooded, and a strait between Denmark and southern Sweden, The Sound, separated these two countries.
The colonization pattern of flora and fauna after LGM was mainly influenced by changes linked to these climatic changes during the late Quaternary (Hewitt, 2004), and the expansions and contractions in the distribution of flora and fauna caused demographic, and spatial events that influenced the pattern of genetic diversity (Hewitt, 2004).
The post-LGM expansion in Europe described by Hewitt (1999) started from southern and eastern refugia with plants and animals moving northwards following a series of rapid temperature increases. However, Pedreschi et al. (2019) suggested that individual species reacted more niche-specific to climate oscillations, thus questioning the generalization of this pattern across species by implying the existence of additional refugia further north or west in Europe depending on the species in question, that is Microtus arvalis (Pedreschi et al., 2019) and Clethrionomys glareolus (Marková et al., 2020). Nevertheless, in most species, genetic diversity is generally higher in populations in the southern areas, as compared with the more recently colonized northern regions. This pattern follows the central-marginal hypothesis predicting a decrease in genetic F I G U R E 1 (a) North-western Europe about 16,000 years ago (above) and North-western Europe about 14,000 years ago (below; Aaris-Sørensen, 2007). (b) the sampling locations of the Northern Birch Mouse analyzed in the study diversity with distance from the glacial refugium (Eckert et al., 2008) and is supported by many phylogeographical studies of animal and plant species.
Until recently, phylogeographic studies primarily relied on the analysis of short fragments or single genes from the mitochondrial genome (mitogenome). Analyzing the full mitogenome (ca. 16-18,000 bp long (Ballard & Whitlock, 2004)) increases the genetic resolution, which is especially important when analyzing closely related populations where diagnostic variation may be limited (Jacobsen et al., 2012;Morin et al., 2010). Most studies have assumed that variation evolves through neutral processes, as the mitogenome is prone to genetic drift given its small effective size (N e ; approximately four times smaller than its nuclear counterparts), due to a haploid nature and maternal inheritance (Ballard & Whitlock, 2004;Birky et al., 1983;Galtier et al., 2009). Genetic drift in combination with strong purifying selection are the two predominant evolutionary forces within the mitogenome. Nevertheless, positive selection might still act on individual genes or single codon positions (Ballard & Whitlock, 2004;Galtier et al., 2009). Several studies have shown evidence for positive selection in a wide array of taxa such as fish and mammals (e.g., Belanger-Deschenes et al., 2013;da Fonseca et al., 2008;Fontanillas et al., 2005;Gagnaire et al., 2012;Garvin et al., 2011;Jacobsen et al., 2014Jacobsen et al., , 2015. Such changes may be linked to metabolic differences and ultimately cover adaptations that are vital for the survival of the specific populations or species. Positive selection may also influence divergence between populations and hence impede reliable estimation of divergence time. This study focuses on the phylogeographic pattern and colonization history of the Northern Birch Mouse (Sicista betulina), (Dipodidae: Zapodinae). This species was among the first rodent immigrants to arrive in Northern Europe and subfossil remains from Denmark date to the Late Glacial period, ~14,000 years ago, when the land was covered by a tundra-boreal forest-like vegetation (Aaris-Sørensen, 1995). The Sicista genus shows highly variable karyotypes and cluster into five major species groups: betulina, caucasica, caudate, tianschanica, and concolor (Lebedev, Poplavskaya, et al., 2020). The betulina group is divided into a Northern species group consisting of S. betulina and Strand's Birch Mouse (S. strandi) and a Steppe group consisting of amongst others, the Southern Birch Mouse (S. subtilis; Holden et al., 2017;Lebedev, Rusin, et al., 2019 Kovalskaya et al., 2011;Lebedev, Rusin, et al., 2019).
The origin of the genus Sicista can be referred to Asia (Kimura, 2011), but in Europe, fossil remains in southeastern France suggest that S. betulina probably appeared during Middle Pleistocene, around 781,000 to 126,000 years ago (Rofes et al., 2012). Further fossil finds of S. betulina are reported from Jersey-La Cotte de St.
Brelade from around Middle-to Late Pleistocene in the Saalian and Weichselian layers (380,000-11,500 years ago; Chaline & Brochet, 1986) and from Boxgrove, West Sussex dating ~Middle Pleistocene (around 550,000-300,000 years ago; AHOB database -Pettitt & White, 2012). These findings suggest that the Northern Birch Mouse may have been present in Northern Europe during interglacial periods, with its range being retracted during the Weichselian glaciation to refugia further south and east within the tundra-forest habitats. Fossil records are found in northern Denmark-in the Nørre Lyngby excavation-that can be dated to ~11,500 BC (Aaris- Sørensen, 1995). In Germany and Hungary, fossil records are dated to ~Middle Pleistocene (around 781,000 to 126,000 years ago) supporting the finds from Jersey and West Sussex, and in Austria, Poland, Switzerland, Romania, Italy, and Spain, the oldest fossils are recorded from Late Pleistocene (126,000-11,430 years ago;Rofes et al., 2012). Given the presumably long persistence of Northern Birch Mouse in Europe, it is likely that several refugia existed during the latest ice age, as observed for many other mammals (Hewitt, 1999(Hewitt, , 2004.
The climatic oscillations following the LGM which led to several glacial comebacks might have prompted several colonization waves of the Northern Birch Mouse into Northern Europe, including Denmark. These may have represented expansions from different refugia with different evolutionary histories, which might explain the observed highly fragmented distribution of the species today.
If so, this is expected to be reflected in the genetic and phylogeographic patterns of the populations.
Here, we investigated the phylogeographical patterns of the Northern Birch Mouse in Northern Europe with special focus on Denmark. To do this, we analyzed full mitochondrial genome sequences (mitogenomes), in combination with genetic variation found in two nuclear genes: beta-myosin heavy-chain (MYH6) and the entire intron 7 of β-fibrinogen (FGB). The samples covered two areas in Denmark (north and south), Norway, Sweden, Estonia, Latvia, Russia, and Slovakia and included an individual of Strands Birch Mouse (S. strandi), a sister species to the Northern Birch Mouse (Cserkész et al., 2015). The colonization history was investigated analyzing the phylogenetic relationship between the mitogenomes and nuclear haplotypes amongst the different sampling locations, estimating divergence time between major phylogeographical clades and populations. Selection was further analyzed to test for deviations from neutrality, and to assess whether adaptive variation existed between populations and species. For the Danish Birch Mouse, which are separated into a northern and southern population (Jensen & Møller, 2007), we tested whether the present populations represented one or more refugia by estimating divergence time using IMa software (Hey & Nielsen, 2007), and examining the population structure. Specifically, we ). This was conducted by an iterative process, that is Sangersequencing and re-designing primers. To avoid numts, the complete mitochondrial genome was amplified using three long fragments of more than ~3 kb (Richly & Leister, 2004). The PCRs were conducted in a total volume of 20 µl using GOTaq Long PCR master mix (Promega) using the thermal profile: 95°C for 15 min, 40 cycles of 94°C for 45 s, ann. 48°C (primer3+6) for 45 s, extension of 72°C for 7 min, and finally extension of 72°C for 10 min. Annealing temperature for primer 2 + 9 was 55°C and identical extension time. For primer 7 + 8, annealing temperature was 66°C and extension time was 72°C for 3 min. The nuclear genes were amplified according to Matocq et al. (2007) and Sanger-sequenced at MACROGEN, Europe.
PCR fragments for mitogenomic analysis were cleaned and concentrations measured using a Qubit 2.0 Flourometer. Libraries were subsequently prepared using NEBNext ® Fast DNA Fragmentation & Library Prep Set for Ion Torrent following the manufacturer's recommendations. Libraries were tagged using Ion express barcodes (1-24) for multiplexing and sequencing was conducted using Ion express template 200 chemistry (paired end technology) on 314 and 316 chips according to the manufacturer's protocol. After assembly, not all mitogenome sequences were complete; hence, primers were designed using the mitogenomes as template filling the gaps using Sanger-sequencing (MACROGENE Europe).

| Mitogenome assembly and annotation
Sequences were assembled using Sequencher ® 5.2.3 (Gene Codes Cooperation) and complete mitochondrion sequence from Microtus rossiaemeridionalis as baseline sequence to help de novo assembly of the sequences in the first step. Reads with a coverage ranging from 10 to 300 were used as the total number of reads from each individual varied. The assembly was conducted choosing the long-

| Nuclear genes and haplotype estimation
To identify the haplotypes of the two nuclear genes, MYH6 and FGB, a Bayesian statistical method, was used implemented in PHASE (Scheet & Stephens, 2006;Stephens et al., 2001) and in DnaSP v 5 (Librado & Rozas, 2009). The identified haplotypes were based on three independent runs for each gene using 1000 iterations.

| Detection of selection
Positive selection acting on the 13 coding mitochondrial genes and the two nuclear genes was tested using the MEME (Mixed Effects Model of Evolution; Murrell et al., 2012) and FUBAR (Fast Unbiased Bayesian AppRoximation; Murrell et al., 2013). These tests are codon-based selection tests implemented in HyPHY (DataMonkey server, (https://www.datam onkey.org)). The methods test for episodic selection (MEME) or diversifying (positive) and negative selection (FUBAR) acting on individual codons. Default settings were used. Only codons shared among ≥2 individuals were used to ensure that potential sequencing errors did not affect the results, and their positions were visualized in the consensus phylogenetic tree estimated in *BEAST (see below).

| Population demography, haplotype networks
Historical demographic population fluctuations were explored using the neutrality of mutations in the five different datasets, tested using Tajima's D test of selective neutrality (Tajima, 1989) and Fu's F S (Fu, 1997). Excess numbers of low frequency  (Bandelt et al., 1999). The network was generated using DnaSP (Librado and Rojas 2009) and POPART with ε set at zero (Leigh & Bryant, 2015).

| Mutation rates estimation
Due to the ongoing debate concerning the use of mutation rates estimated from evolutionary old time points in the speciation history and the applicability of these rates to estimate divergence times between populations within a species (Ho et al., 2005;Subramanian & Lambert, 2011), two substitution-rates were used to date observed splitting events. Consequently, divergence time was estimated based on a fast 7.3%/site/MYR , and a slow 2%/site/MYR (Cheng et al., 2019) CytB mutation rates for Dipodidae. Given these two mutation rates were estimated for CytB, we used the formula by Andersen et al. (2017) to estimate the mutation rates for the full mitogenome. The calculations were based on theta (θ) estimated in DNAsp 6.10.03 (Rozas et al., 2017), the assumption of no recombination and identical effective female population size (N fe = θ/2μ) between the regions: μgene = ((θgene × μmitogenome)/θmitogenome)). The new estimated mutation rate for the full mitogenome based on the fast mutation rate was 5.813%/site/Myr, and 1.519%/site/Myr for the slow mutation rate.

| Population structure
Population structure was analyzed using a pairwise Ф ST test conducted on pairwise nucleotide differences in ARLEQUIN v3.5.2 (Excoffier & Lischer, 2010). The tests were run for 10,000 permutations over individual haplotypes among populations. Furthermore, we applied the non-parametric method (Kolmogorov-Smirnov, (K-S) test) suggested by Carr et al. (2015) that analyze "the cumulative pairwise mismatch distribution between populations." (Carr et al., 2015) The Kolmogorov-Smirnov test does not assume that data are sampled from a specific distribution, testing the null hypothesis that data originate from the same population with identical distribution. This was performed to analyze the phylogenetic and phylogeographic structure among the regions.
The nucleotide variants from all mitogenomes within a geographical region (i.e., DK_North) were compared pairwise to the nucleotide substitution sites from another region (i.e., DK_South) based on the aligned mitogenomes. The cumulative frequency curve was formed for each and the largest, unspecified observed difference, D (the vertical offset of the curves), over the range of all differences between the cumulative frequency distributions was identified. The significance of the obtained maximum, D, of the pairwise comparison was assessed by Da, the critical value which has been shown (Kolmogorov, 1933) to be independent of the reference distribution. All the pairwise Kolmogorov-Smirnov tests were performed using the Web-based statistical software, http://www.physi cs.csbsju.edu/stats/ KS-test.html. Due to very low sample size in three of the geographical areas (N < 5, Russia, Norway and Latvia), this analysis was only applied to the remaining six regions.

| Phylogenetic analysis and divergence time estimation
Phylogenetic analysis was conducted on the mitogenome dataset to explore the evolutionary relationships and estimate Time of Most Recent Common Ancestor (TMRCA). First, to investigate the overall evolutionary relationships, a Neighbour joining (NJ) phylogeny was built in MEGA v 6.06 (Tamura et al., 2013) using the TrN substitution model and 1000 bootstraps. Then, Bayesian phylogenetic analysis and estimation of TMRCA was performed using BEAST v1.8.4 (Drummond & Rambaut, 2007). A concatenated dataset including all 13 mitochondrial coding genes was analyzed. The sequences were further divided into the three codon positions to allow separate estimation of mutation rate. The best substitution model was chosen based on a maximum likelihood search conducted in MEGA. Initial analysis found the HKY + G + I model to constitute the most accurate substitution model for the dataset. The final analysis was based on the simpler HKY model due to problems with MCMC convergence in BEAST when using the HKY + G + I model. However, the phylogeny and the estimated TMRCA of all major groups were extremely similar between analyses supporting a limited effect of the substitution models. Initial phylogenetic analysis showed the presence of three major clades. These may have different demographic histories, which could affect divergence estimates if not accounted for (Ho et al., 2008). Thus, to allow differences in population size changes, the final analysis was conducted in *BEAST (Heled & Drummond, 2010). This program uses two different tree priors: a speciation prior for between-clade branch patterns and a coalescence prior for the within-clade branch patterns, which allows inferences of different demographic history between lineages.
Piecewise linear and constant root, piecewise linear and piecewise constant change models of effective population size were used as coalescence priors to compare the robustness of the divergence estimates given differences in demographic history.
The YULE prior was used as the between-clade prior. To estimate the best set of models, bayes factor calculation, as proposed by Suchard et al. (2001), was conducted in BEAST. All phylogenetic analyses used a constant clock as a preliminary test in BEAST showed no significant rate heterogeneity among branches. Both the slow and the fast mutation rates were used to estimate divergence time in years.
The final MCMC samples were based on a run for 50,000,000 generations, and genealogies were sampled every 5000 generations with 10% discarded as burn-in. Examination of convergence and effective sample size (ESS) values was conducted using TRACER v1.5 (Drummond & Rambaut, 2007). All parameters had ESS values >200, and additional runs gave similar results. The maximum clade credibility tree with mean heights for branches was estimated in the program TREEANNOTATOR (Drummond & Rambaut, 2007) with 10% burn-in and visualized and edited in the program FIGTREE v1.3.1 (Andrew Rambaut, University of Edinburgh, http://tree.bio.ed.ac. uk/softw are/figtr ee/). Divergence time between the different geographical groups was estimated using IMa software (Hey & Nielsen, 2007) using the full mitogenomes. Only groups with more than five individuals were analyzed to avoid spurious results due to low sample size. IMa infer divergence time between two populations from a common ancestral population. It assumes constant population sizes and no intragenic recombination, which is not occurring in mtDNA. A model with isolation without migration (m=0) was chosen to mimic founding of the Danish populations. Population sizes (q) were set to 300 and divergence time, t, to 20. A geometric heating scheme with 15 MCMC chains with the parameters, g1 = 0.7, g2 = 0.9, L = 3.0 and a burn-in = 30,000 was used to explore the parameter space.
Convergence was evaluated by ESS and performing two runs for each pairwise estimation of divergence time showing similar results.
The dataset was analyzed using the HKY substitution model and an inheritance scalar of 0.25 to account for the lower effective size of mtDNA compared with nuclear DNA. As in *BEAST, the fast and slow mutation rates were used for divergence time estimation in years.

| Haplotype variation
Fifty-four mitogenomes, 69 FGB and 66 MYH6 nuclear sequences of S. betulina and one S. strandi were used in the analysis (Table 1). In total, the 54 different mitogenome haplotypes were detected representing 1699 segregating sites. For the two nuclear genes, seven different sequences were observed in MYH6 with six segregating sites, and in FGB, ten sequences and eleven segregation sites were found after specifying the haplotypes in PHASE including S. strandi in both analyses (Scheet & Stephens, 2006;Stephens et al., 2001).
The genetic diversity in terms of nucleotide diversity (π) among the mitogenome samples ranged from ~0.614% ± 0.144% observed in the DK_South (DKS) sample to ~3.728% ± 0.783% observed in the Estonia sample (Table 1). For the nuclear genes, nucleotide diversity for MYH6 ranged from ~0.316% ± 0.038% in DK_North (DKNT) to ~0.553% ± 0.169% in the Latvian sample. In the FGB gene, the nucleotide diversity ranged from ~0.053% ± 0.021% in DK_North and to ~0.282% ± 0.102% in the sample from Latvia. Sample size was, however, very low from Latvia.  (Grantham, 1974). No significant selection was detected acting in the two nuclear genes (data not shown).

| Population demography, haplotype relationship
Exploring the possibility of detecting signs of population expansion or bottlenecks, Tajima's D and Fu's F S were estimated for all datasets ( Table 1). None of the two parameters were significant in any of the datasets used. This supported the assumption of constant population sizes, which is an essential assumption for the divergence time In FGB, 10 haplotypes were detected. SBdkFGB1 was found in all locations while SBdkFGB2 was detected in all locations except Norway, Estonia, and Slovakia. SBdkFGB5 and SBdkFGB6 were TA B L E 1 Sample size (N), Hd (haplotype diversity), π (nucleotide diversity), SD (standard deviation), Tajima D and Fu's F s (Fu, 1997)

| Population structure
Despite low sample sizes, the pairwise analysis of molecular variance  (Table 3; Figure 2).

| Phylogeny and divergence time estimation
The three different Bayesian analyses showed very similar estimates of TMRCA. However, the analysis using the piecewise linear and constant root coalescence prior showed the best fit to the data (bayes factor >9000) and were therefore used for visualization of evolutionary relationships and reporting of TMRCA ( Figure 3)  BC-5314 BC) using the fast mutation rate ( Table 4; Figure 4). The different estimates should, however, be interpreted cautiously given the uncertainty regarding low sample sizes and mutation rates. This is furthermore reflected in the width of the peaks in the probability of divergence ( Figure 4) and the 90% high point density interval (HPD , Table 4).

| DISCUSS ION
This study is the first to analyze the phylogeographic history of the

| Genetic variation
Genetic diversity was higher for the mitogenomes compared with the nuclear markers. This is expected as mitochondrial DNA in general evolves at a faster rate than nuclear genes and hence accumu- Note: Bold = significant after sequential Bonferroni correction, α = 5% (Rice, 1989).

F I G U R E 2
Cumulative pairwise distribution curves for the Kolmogorov-Smirnov test based on the nucleotide differences between the mitogenomes sampled in the seven locations of Northern Birch Mouse. Norway and Russia are included for illustration despite the low sample sizes during Paleocene/Early Eocene or might be influenced by introgression from other closely related species (Marivaux et al., 2004;Pisano et al., 2015;Yue et al., 2015).
The central-marginal hypothesis predicts a decrease in genetic diversity with increasing distance from the glacial refugium (Eckert et al., 2008). This implies that the Danish and Swedish locations are expected to hold a lower nucleotide diversity compared with more southern and eastern locations such as Slovenia and Estonia. Despite the low sample sizes, this tendency was tentatively suggested by the observed variation in both the mitogenomes and the nuclear genes.
As such, this study suggests recolonization of northern Europe from more southern (Carpathian) or eastern areas that served as refugia during the latest ice age. This is supported by many phylogeographical studies of steppe fauna which find higher genetic diversity in Southern and Eastern Europe (Kajtoch et al., 2016;Marková et al., 2020).  TA B L E 4 IMA results estimated without migration (Hey & Nielsen, 2007) shown for the samples with more than five individuals

| Selection and mutation rate
The 13 mitochondrial coding genes were observed mainly to be under purifying selection as observed in other studies (Andersen et al., 2017;Jacobsen et al., 2015;Tomasco & Lessa, 2011). However, some sites showed significant dN/dS-values >1, suggesting possible positive selection. These non-synonymous sites were mainly found in the ND genes, which in general show higher dN/dS than other mitochondrial genes (Andersen et al., 2017), and almost exclusively included situations where amino acid substitutions occurred independently in multiple places in the phylogenetic tree. ND genes belong to the electron transport system I (ETS I) in the four components of the mitochondrial ETS and have been observed having a faster mutation rate compared with the other three ETS (Matosiuk et al., 2014;Zhang & Broughton, 2013). This may be the result of functional constraints, as ND genes probably are of less importance for metabolic output than the other mitochondrial genes, which may lead to relaxed purifying selection (Blier et al., 2006). Positive selection was not further supported as the amino acid replacements showed very similar physiochemical properties, which suggests either neutral or near neutral changes. The same pattern was also observed in the four other genes with detected candidate codon positions for positive selection (COX1, COX3, CytB, and ATP6). Hence, while these could cover positive selected sites, we find it more parsimonious that the observations are a result of relaxed purifying selection and reflect neutrally or slightly deleterious mutations.
Evidence for this has been found in many other mitogenome studies and may relate to change in the efficiency of purifying selection (Ho et al., 2005;Jacobsen et al., 2015), which might need time to remove slightly deleterious mutations. As a result, the observed mutation rate may also be higher, which supports the use of the fast mutation rate for divergence estimation. The rate was estimated to 5.81% substitutions/site/Myr based on a CytB mutation rate of 7.3%/site/Myr in rates for Dipodidae . This rate is slightly faster than the mitogenome mutation rate of 3.2%/ site/Myr observed in M. musculus domesticus from the Kerguelen Archipelago (Hardouin & Tautz, 2013).

| Population structure and demographic history
The prevalent phylogenetic methods used in population genetic studies have problems when all individuals are genetically different as in this study. The pairwise AMOVA based on differences between haplotypes/mitogenomes (Ф ST ) among the populations is too conservative as all the variance is observed among the individuals.
Consequently, the population structure pattern is difficult to resolve and only DKS was found significantly different from the more distant localities (Estonia, Slovakia, Russia). Using nuclear markers with slower mutation rates on the other hand revealed too little variation. Only FBG was able to detect significant genetic differentiation between the more distant locations, Slovakia and Estonia, and the Danish locations (except Slovakia-DKS). Due to the very high mitogenome variation, we chose to use the non-parametric K-S test to analyze the population relationship and phylogenetic history as suggested by Carr et al. (2015). These authors used the

| Phylogeograpic patterns and divergence times
In general, the genetic relationship reflected by the median-joining haplotype network based on mitogenomes supported a closer relationship between individuals within the different locations compared to between locations. However, the huge number of mutations separating the individuals even within the same geographical area support the hypothesis that colonization of Northern Europe might have occurred from different refugia. This combined with limited mobility of the species, repeated changes in climate and vegetation together with rapid radiation connected to Rodentia are probably the drivers behind the observed genetic pattern.
The nuclear haplotype network reflected a very close relationship between all locations illustrating the limited observed variation in those genes prompt by an expected lower nuclear mutation rate and shorter sequences leading to a lower genetic resolution.
Nevertheless, both genes showed patterns in support of expansion from several refugia.
The phylogenetic tree estimated TMRCA between the S. strandi clade and the Northern Birch Mouse clades to ~527,033 years ago (494,207-561,346 95%HPD) using the fast substitution rate (µ fast ).
This is slightly more recent compared with the divergence time, ~670,000 years, suggested by Rusin et al. (2018) based on variation in CytB and five nuclear genes but earlier than the around 1 Mya (Early Pleistocene) estimated by Lebedev, Rusin, et al. (2019) and well within the suggested appearance of S. betulina from the fossil record in southeastern France ~781,000 to 126,000 years ago (Rofes et al., 2012). The differences in the TMRCA estimates can probably be attributed to the different kind and number of markers (nuclear and/or mitochondrial genes), and mutation rates used in the different studies. Common for these studies are that they support the species status of S. betulina and S. strandi. Interestingly, we found that three individuals from Slovakia and Estonia actually Swedish populations (Bennike et al., 2012).

Divergence estimation between DK_North and DK_South
suggested divergence during Early Neolithic to Late Iron Age (400 AD-1050 AD) depending on mutation rate (Odgaard & Rasmussen, 2000). In general, the landscape was influenced by anthropogenic activities and climatic changes during these time periods causing changes in forest land cover, the type of forest and the surrounding land use, increasing the grassland and heathland cover which probably have determined the distribution of the Northern Birch Mouse (Nielsen & Odgaard, 2010). However, the result remains tentative as the small sample size prevents accurate analysis of population structure and divergence timing. Thus, whether the two Danish areas indeed represent two recently reproductively isolated populations cannot be fully demonstrated by this study and we therefore recommend additional genetic studies using larger sample sizes and potentially more genetic markers to investigate this possibility.

ACK N OWLED G M ENTS
A special thanks to all the researchers contributing with material to

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.